{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualize ICESat-2 ATL07 Data\n",
    "\n",
    "ICESat-2 hackweek final project  \n",
    "June 16, 2020  \n",
    "Nicole Abib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import datetime as dt\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import h5py\n",
    "from scipy.interpolate import interp1d\n",
    "from astropy.time import Time\n",
    "import sys  \n",
    "sys.path.insert(0, '/home/jovyan/leading_to_phytoplankton/scripts')\n",
    "import readers as rd\n",
    "\n",
    "# magic function to enable interactive plotting\n",
    "%matplotlib widget "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    }
   ],
   "source": [
    "# Import ATL07 Data\n",
    "beam='gt1l' # Beam of interest\n",
    "data_loc='/home/jovyan/shared/leading_to_phytoplankton/'\n",
    "is2_fn = 'IS2_S2/ATL07-01_20190805213851_05930401_002_02.h5'\n",
    "is2_f = h5py.File(data_loc+is2_fn, 'r')\n",
    "df07= rd.getATL07(is2_f,beam)\n",
    "df07.head()\n",
    "print(is2_f['orbit_info/sc_orient'][0]) # Check forward/backward orientation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lats</th>\n",
       "      <th>lons</th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "      <th>heights</th>\n",
       "      <th>dt</th>\n",
       "      <th>conf</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>79.993666</td>\n",
       "      <td>40.449280</td>\n",
       "      <td>8938354.0</td>\n",
       "      <td>3306.422607</td>\n",
       "      <td>-2426.635986</td>\n",
       "      <td>5.027759e+07</td>\n",
       "      <td>-2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>79.993605</td>\n",
       "      <td>40.450261</td>\n",
       "      <td>8938344.0</td>\n",
       "      <td>3289.306152</td>\n",
       "      <td>234.947739</td>\n",
       "      <td>5.027759e+07</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>79.993611</td>\n",
       "      <td>40.450256</td>\n",
       "      <td>8938344.0</td>\n",
       "      <td>3289.260986</td>\n",
       "      <td>242.102203</td>\n",
       "      <td>5.027759e+07</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>79.993612</td>\n",
       "      <td>40.450246</td>\n",
       "      <td>8938344.0</td>\n",
       "      <td>3289.424561</td>\n",
       "      <td>216.635086</td>\n",
       "      <td>5.027759e+07</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>79.993617</td>\n",
       "      <td>40.450247</td>\n",
       "      <td>8938345.0</td>\n",
       "      <td>3289.272705</td>\n",
       "      <td>240.364456</td>\n",
       "      <td>5.027759e+07</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        lats       lons          x            y      heights            dt  \\\n",
       "0  79.993666  40.449280  8938354.0  3306.422607 -2426.635986  5.027759e+07   \n",
       "1  79.993605  40.450261  8938344.0  3289.306152   234.947739  5.027759e+07   \n",
       "2  79.993611  40.450256  8938344.0  3289.260986   242.102203  5.027759e+07   \n",
       "3  79.993612  40.450246  8938344.0  3289.424561   216.635086  5.027759e+07   \n",
       "4  79.993617  40.450247  8938345.0  3289.272705   240.364456  5.027759e+07   \n",
       "\n",
       "   conf  \n",
       "0    -2  \n",
       "1     2  \n",
       "2     2  \n",
       "3     2  \n",
       "4     2  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Import ATL03 Data from the same RGT for comparison\n",
    "is2_fn = 'IS2_S2/ATL03_20190805215948_05930404_002_02.h5'\n",
    "is2_f = h5py.File(data_loc+is2_fn, 'r')\n",
    "df03=rd.getATL03(is2_f, beam)\n",
    "df03.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Look at a smaller section"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "lonmin,lonmax,latmin,latmax= [-180,180,82,87.5]\n",
    "region=(df07.lats>latmin)&(df07.lats<latmax)&(df07.lons>lonmin)&(df07.lons<lonmax)\n",
    "df07_cut=df07[region]\n",
    "# cut df03 to region of interest\n",
    "region=(df03.lats>latmin)&(df03.lats<latmax)&(df03.lons>lonmin)&(df03.lons<lonmax)\n",
    "df03_cut=df03[region]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Apply vertical corrections to ATL03 to allow comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lats</th>\n",
       "      <th>lons</th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "      <th>heights</th>\n",
       "      <th>dt</th>\n",
       "      <th>conf</th>\n",
       "      <th>correction</th>\n",
       "      <th>height_corr</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>815511</th>\n",
       "      <td>82.000008</td>\n",
       "      <td>37.311930</td>\n",
       "      <td>9168910.0</td>\n",
       "      <td>3290.816650</td>\n",
       "      <td>18.217131</td>\n",
       "      <td>5.027762e+07</td>\n",
       "      <td>1</td>\n",
       "      <td>22.380521</td>\n",
       "      <td>-4.163390</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>815512</th>\n",
       "      <td>82.000032</td>\n",
       "      <td>37.311886</td>\n",
       "      <td>9168912.0</td>\n",
       "      <td>3290.770996</td>\n",
       "      <td>26.980247</td>\n",
       "      <td>5.027762e+07</td>\n",
       "      <td>1</td>\n",
       "      <td>22.380583</td>\n",
       "      <td>4.599665</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>815513</th>\n",
       "      <td>82.000039</td>\n",
       "      <td>37.311874</td>\n",
       "      <td>9168913.0</td>\n",
       "      <td>3290.767090</td>\n",
       "      <td>27.991995</td>\n",
       "      <td>5.027762e+07</td>\n",
       "      <td>1</td>\n",
       "      <td>22.380598</td>\n",
       "      <td>5.611397</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>815514</th>\n",
       "      <td>82.000082</td>\n",
       "      <td>37.311786</td>\n",
       "      <td>9168918.0</td>\n",
       "      <td>3290.833008</td>\n",
       "      <td>20.700541</td>\n",
       "      <td>5.027762e+07</td>\n",
       "      <td>1</td>\n",
       "      <td>22.380707</td>\n",
       "      <td>-1.680166</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>815515</th>\n",
       "      <td>82.000100</td>\n",
       "      <td>37.311756</td>\n",
       "      <td>9168920.0</td>\n",
       "      <td>3290.747314</td>\n",
       "      <td>35.260612</td>\n",
       "      <td>5.027762e+07</td>\n",
       "      <td>0</td>\n",
       "      <td>22.380753</td>\n",
       "      <td>12.879859</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             lats       lons          x            y    heights            dt  \\\n",
       "815511  82.000008  37.311930  9168910.0  3290.816650  18.217131  5.027762e+07   \n",
       "815512  82.000032  37.311886  9168912.0  3290.770996  26.980247  5.027762e+07   \n",
       "815513  82.000039  37.311874  9168913.0  3290.767090  27.991995  5.027762e+07   \n",
       "815514  82.000082  37.311786  9168918.0  3290.833008  20.700541  5.027762e+07   \n",
       "815515  82.000100  37.311756  9168920.0  3290.747314  35.260612  5.027762e+07   \n",
       "\n",
       "        conf  correction  height_corr  \n",
       "815511     1   22.380521    -4.163390  \n",
       "815512     1   22.380583     4.599665  \n",
       "815513     1   22.380598     5.611397  \n",
       "815514     1   22.380707    -1.680166  \n",
       "815515     0   22.380753    12.879859  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# mean sea surface correction - interpolating time with corrections\n",
    "mss_corr=interp1d(df07.dt.values,df07.mss.values)\n",
    "# inverted barometer correction\n",
    "ib_corr=interp1d(df07.dt.values,df07.ib.values)\n",
    "# ocean tide correction\n",
    "ocean_corr=interp1d(df07.dt.values,df07.ocean.values)\n",
    "# long period equilibrium tide correction\n",
    "lpe_corr=interp1d(df07.dt.values,df07.lpe.values)\n",
    "\n",
    "# Apply corrections  at sampling rate of dt from atl03\n",
    "df03_cut['correction']=mss_corr(df03_cut.dt.values)+ib_corr(df03_cut.dt.values)+ocean_corr(df03_cut.dt.values)+lpe_corr(df03_cut.dt.values)\n",
    "df03_cut['height_corr']=df03_cut.heights-df03_cut.correction # subtract correction from ATL03 heights\n",
    "df03_cut.head()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "eebbe64cbb60433cb157aa402e54d105",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "var= 'stype' \n",
    "\n",
    "## we will want to set colorbar parameters based on the chosen variable\n",
    "vmin=0\n",
    "vmax=10\n",
    "ticks=np.arange(vmin,vmax,1)\n",
    "\n",
    "plt.figure(figsize=(8,8), dpi= 90)\n",
    "ax = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0)) # choose polar sterographic for projection\n",
    "ax.coastlines(resolution='50m', color='black', linewidth=1)\n",
    "# ax.set_extent([-80,-40,80,86], ccrs.PlateCarree())\n",
    "# ax.set_extent([-180,180,82,87.5], ccrs.PlateCarree())\n",
    "ax.set_extent([-180, 180, 60, 90], ccrs.PlateCarree())\n",
    "plt.scatter(df07_cut['lons'][::100], df07_cut['lats'][::100],c=df07_cut[var][::100], cmap=plt.cm.get_cmap('RdYlGn').reversed(), vmin=vmin,vmax=vmax,transform=ccrs.PlateCarree())\n",
    "plt.colorbar(label=var, shrink=0.5, ticks=ticks,extend='both');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0e85b032d6ef422a92d9b23cd623e5e2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax=plt.subplots(5,1,figsize=(15,20))\n",
    "\n",
    "ax[0].scatter(df07_cut.lats,df07_cut.photon_rate,s=2,c='gray', label='ATL07 all surface type')\n",
    "ax[0].scatter(df07_cut[df07_cut.stype>1].lats,df07_cut[df07_cut.stype>1].photon_rate,s=10,c='b', label= 'ATL07 surface type lead')\n",
    "ax[0].legend(fontsize=15)\n",
    "ax[0].set_ylim(0,20)\n",
    "ax[0].set_ylabel('ATL07 \\nphoton rate',fontsize=15)\n",
    "\n",
    "ax[1].scatter(df07_cut.lats,df07_cut.gauss,s=2,c='gray', label='ATL07 all surface type')\n",
    "ax[1].scatter(df07_cut[df07_cut.stype>1].lats,df07_cut[df07_cut.stype>1].gauss,s=10,c='blue',label= 'ATL07 surface type lead')\n",
    "ax[1].legend(loc=1,fontsize=15)\n",
    "ax[1].set_ylabel('ATL07 \\ngaussian width',fontsize=15)\n",
    "ax[1].set_ylim(0,.6)\n",
    "\n",
    "ax[2].grid()\n",
    "ax[2].scatter(df07_cut.lats,df07_cut.stype,s=2,c='gray',label='ATL07 all surface type')\n",
    "ax[2].scatter(df07_cut[df07_cut.stype>1].lats,df07_cut[df07_cut.stype>1].stype,s=10,c='blue',label= 'ATL07 surface type lead')\n",
    "ax[2].legend(fontsize=15)\n",
    "ax[2].set_ylim(-1,10)\n",
    "ax[2].set_ylabel('ATL07 \\nsurface type',fontsize=15)\n",
    "ax[2].set_yticks(np.arange(0,10,1))\n",
    "\n",
    "ax[3].scatter(df03_cut.lats[df03_cut.conf<3],df03_cut.height_corr[df03_cut.conf<3],c='r',s=.5,alpha=.3)\n",
    "ax[3].scatter(df03_cut.lats[df03_cut.conf>2],df03_cut.height_corr[df03_cut.conf>2],c='g',s=.5,alpha=.3)\n",
    "ax[3].scatter(0,0,c='g',s=5,alpha=1, label = 'ATL03 med/high conf')\n",
    "ax[3].scatter(0,0,c='r',s=5,alpha=1, label= \"ATL03 low conf\")\n",
    "\n",
    "ax[3].scatter(df07_cut.lats,df07_cut.heights,marker='s',s=5,c='k',alpha=.5, label= 'ATL07 heights')\n",
    "ax[3].scatter(df07_cut[df07_cut.stype>1].lats,df07_cut[df07_cut.stype>1].heights,marker='s',s=10,c='b', label= 'ATL07 surface type lead')\n",
    "ax[3].legend(ncol=2, loc=4,fontsize=15)\n",
    "ax[3].set_ylim(-1,2)\n",
    "ax[3].set_ylabel('ATL03 and ATL07 \\nsurface heights',fontsize=15)\n",
    "\n",
    "ax[4].scatter(df07_cut.lats,df07_cut.ssh_flag,c='k',s=2,marker='.',label='ATL07 all surface type')\n",
    "ax[4].scatter(df07_cut[df07_cut.stype>1].lats,df07_cut[df07_cut.stype>1].ssh_flag,c='b',s=10,label='ATL07 surface type lead')\n",
    "ax[4].legend(fontsize=15)\n",
    "ax[4].grid()\n",
    "ax[4].set_ylim(-1,2)\n",
    "ax[4].set_yticks(np.arange(0,2,1))\n",
    "ax[4].set_yticklabels(['sea ice', 'sea surface'],fontsize=12)\n",
    "ax[4].set_ylabel('ATL07 ssh flag',fontsize=15)\n",
    "\n",
    "ax[4].set_xlabel('latitude',fontsize=15)\n",
    "\n",
    "for a in np.arange(0,5):\n",
    "    ax[a].set_xlim(latmin,latmax)\n",
    "#     ax[a].axvline(df07_cut[df07_cut.ssh_flag==1].lats.values[0],c='y',alpha=.2, linewidth=10)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
